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Abstract 

We suggest a numerical integration procedure for solving the equations of motion of 
certain classical spin systems which preserves the underlying symplectic structure 
of the phase space. Such symplectic integrators have been successfully utilized for 
other Hamiltonian systems, e. g. for molecular dynamics or non-linear wave equa- 
tions. Our procedure rests on a decomposition of the spin Hamiltonian into a sum 
of two completely integrable Hamiltonians and on the corresponding Lie- Trotter 
decomposition of the time evolution operator. In order to make this method widely 
applicable we provide a large class of integrable spin systems whose time evolution 
consists of a sequence of rotations about fixed axes. We test the proposed symplec- 
tic integrator for small spin systems, including the model of a recently synthesized 
magnetic molecule, and compare the results for variants of different order. 
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1 Introduction 



To calculate the time evolution of classical spin systems is an important task in 
condensed matter physics. For example, the cross section of neutron scattering 
at a spin system is proportional to the Fourier transform of the time- depending 
auto-correlation function, see [1], which can often be calculated in the classi- 
cal limit. Completely integrable spin systems are rare, that is, in most cases 
an analytical calculation of the time evolution is not possible and one is lead 
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to employ numerical integration methods. Since classical spin systems arc in- 
stances of Hamiltonian systems, it is advisable to use numerical integrators 
which preserve the underlying symplectic structure of the phase space. Such 
"symplectic integrators" have been considered in the last decades [2] and have 
been applied to a variety of problems, ranging from molecular dynamics [3] to 
the nonlinear Schrodingcr equation [4,5]. 

Unfortunately, symplectic integrators for spin systems have only rarely been 
considered in the literature, see [6,7,8]. The method of the independent time 
evolution of sublattices, proposed in [6,9,10], is volume-preserving but not 
symplectic, see 2.2.1 and [6]. Inspired by [10], we suggest to construct sym- 
plectic integrators based on a splitting of the spin Hamiltonian into two com- 
pletely integrable Hamiltonians belonging to a special kind of systems [11,12]. 
These systems are called "i3-partitioned systems" and their time evolution can 
be calculated as a sequence of rotations about fixed axes [12]. This general- 
izes the Stormer/Verlet scheme based on separable Hamiltonians of the form 

//-r(p) + y(q). 

In section 2 we provide the general definitions and results we need from an- 
alytical mechanics (section 2.1) and from the field of symplectic integrators 
based on Lie- Trotter decompositions of the time evolution operator (section 
2.2). The reader who is not famihar with the differential geometric background 
may skip the technical details and only draw the moral that a symplectic inte- 
grator approximates the exact time evolution by a sequence of calculable time 
evolutions corresponding to auxiliary Hamiltonians. In section 3 we shortly re- 
capitulate the theory of i3-partitioned systems from [12]. In order to test our 
suggestions we have implemented various variants of symplectic integrators 
and applied them to selected small spin systems, see section 4. We report the 
fluctuation of the total energy about its initial value as opposed to the con- 
stant drift for a non-symplcctic Runge-Kutta method (RK4), sec section 4.1. 
For two integrable spin systems we compare the errors of the various symplec- 
tic methods, including RK4, see sections 4.2.1 and 4.2.2. Finally, we compare 
the errors of five symplectic integrators for the integrable N — 5 spin pyramid 
and fixed runtime, see section 4.3. We close with a summary and outlook. 



2 Definitions and general results 



We will only formulate the pertinent definitions for symplectic integrators in 
the context of spin systems. For the general case there are excellent sources 
available in the literature, see e. g. [13,14] for analytical mechanics and [2] for 
symplectic integrators. 
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2. 1 Generalities 



Classical spin configurations can be represented by A'"-tuples of unit 3- vectors 
s = (si, . . . , Sat), = 1 for = 1, . . . , N . The compact, 2A'"-dimensional 

manifold of all such configurations is the phase space of the spin system 



= Piv = {(si, . . . , sjv) I = 1 for /X = 1, . . . , 7V} 



A special coordinate system is given by the 2N local functions 
Vix-i Zfx'-V ^"^1 implicitly defined by 



1 - zlcosipi. 



1 - zlsm^^ 



, /X = 1, . . . , AT . 



(2) 



A tangent vector of P at a point s G P can be represented by an A"-tuple 
t = (ti, . . . ,tiv) of 3- vectors satisfying the constraint 



(3) 



If a, b are two tangent vectors at s e 7^, the assignment 



N 



(4) 



defines a non-degenerate, closed 2-form, that is, a symplectic form u. In the 
coordinate system (2) uj can locally be written in the form 



N 



(5) 



hence {ip^, -2^)^=i,...,Ar are canonical coordinates w. r. t. cu. The volume form 
dV is defined by 



dv = uj^ = uj ^uJ ^uJ 



(6) 



and has the local coordinate representation 
dV — d(pi A dzi A ... A d(fiN A dz^ ■ 



(7) 
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A smooth Hamiltonian : — > R generates the Hamiltonian vector field Xh 
imphcitly defined by 

ix^ uj = u{Xh, )^dH . (8) 



The corresponding Hamiltonian equations of motion are 

I s(t) = Xnm) (9) 



and assume their usual form 

d , , dB. d , , dB. ,^ , , 



in the canonical coordinate system (2). By writing the solution s(t) of (10) in 
the form s(i) = ^t(if)(s(0)) we obtain the Hamiltonian flow Tt{H) -.V^V. 
It is defined for all initial values s(0) and for all t G M since V is compact, 
i. e. Xh is a complete vector field. Analogously, the fiow of a general vector 
field can be defined. 



A smooth map : V ^ V is called symplectic iff it preserves the symplec- 
tic form, i. e. iff (f)*uj = uo. Every symplectic map preserves the phase space 
volume, but not conversely, see the counter-example below. Any Hamiltonian 
flow J^t{H) is symplectic, cf. , for example, theorem 8.1.9 in [15]. Conversely, 
if the flow J^t of a complete vector fleld X is symplectic, then 

£xa;(s) = -^(.F;a;)(s)U = 0, (11) 



where jCx is the Lie derivative. Hence = jCxi^ — ix dw + dixoj = + dixoj , 
that is, ixOJ is a closed 1-form, and has, by the Poincare lemma, locally the 
form ixOJ = dK. To summarize: symplectic flows are, at least locally, generated 
by suitable Hamiltonians K. 



2.2 Symplectic integrators 



From an abstract point of view, a symplectic integrator is an approximation 
of some exact flow J^t{H) by the composition of symplectic maps (p^ : V ^ V, 
which can be calculated analytically or numerically exact. In this article we 
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Table 1 

Various decompositions of the form (16) which give rise to different symplectic 
integrators. 



Name 


Abbr. 


Order 


Coefficients 


Ref. 


Suzuki- Trotter 


STl 


1 


ai = b^ = h 


[161 


Suzuki-Trotter 


ST2 


2 


oi = 02 = 5, 6i = 1, 62 = 


[17] 


Suzuki-Trotter 


ST4 


4 


01=06 = 1, 66 = 

hi = n.o = ho = Ha = n.ci = h^ = Tt 

'-'i "-^z '-'Z '-'4 ^0 IT 

03 = 04 = 63 = 1 - 4p 

^' ~ 4-41/3 


[18] 


Forest-Ruth 


FR 


4 


Oi = O4 = |, 0,2 = 03 = 

6^ = 63 = 6*, 64 = 
^^2 = 1-2^,^ = ^ 


[19] 


Optimized 
Forest-Ruth 


OFR 


4 


ai = 05 = ^, 02 = 04 = X, 65 = 
03 = 1 - 2(^ + x), 61 = 64 = ^ 
62 = 63 = A = -0.09156203 
^ = 0.17208656, x = -0.16162176 


[20] 



assume that the Hamiltonian H is decomposable into completely integrable 
Hamiltonians in the form 

H^Y.Hi (12) 



and that the (pi, are the Hamiltonian flows corresponding to certain Hj. The 
precise form of the correspondence is given by a Lie- Trotter decomposition of 
the flow J-t{H) written as an exponential operator 

J^t{H) = e*^ . (13) 



In order to make sense of (13) we have to linearize the Hamiltonian equations 
of motion. To this end we consider J-t{H) acting on functions f -.V ^ C via 

.F,(i/)7(s) = /(-F_,(if)(s)) . (14) 



If / runs through C'^{V,dV), the Hilbert space of (equivalence classes of) 
square-integrable complex functions, (14) deflnes a continuous, unitary 1- 
parameter group, see section 7.4 of [15] for details. By Stone's theorem, this 
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group has the form (13) with an anti-selfadjoint operator Ti. One can show 
that H. can be expressed by means of the Poisson bracket according to 

Hf = {H, /} = u;{Xh, Xf), f smooth , (15) 

but we will not need this in the sequel. (13) is only needed to provide a basis 
for using the techniques of Lie-Trotter decomposition for Hamiltonian flows. 

For sake of simplicity let us consider the special case H = Hi + H2 and hence 
7i = Hi + 7i2- We are looking for £-th order Lie- Trotter decompositions which 
have the form 

k 

^t{m+m) ^ -Q ^a,tm^b,tm + o{t^+^) . (16) 

1=1 

Both sides of (16) are expanded into power series in terms of t and set equal 
up to terms including t^. This yields a system of, in general, non-linear equa- 
tions for the unknown coefficients Oj, bi. Except for ^ = 1 the corresponding 
solutions are not unique. Hence there exist several decompositions and thus 
several symplectic integrators of the same order £. In this article we will use 
the decompositions enumerated in table 1. All corresponding integrators are 
symmetric, or time-reversible, see [2]. Obviously, the Lie- Trotter decomposi- 
tion (16) is a good approximation only for small t. Therefore the given time 
interval [0, t\ is usually spht into L intervals of length A and (16) is separately 
applied to each time step A. Hence, apart from the choice of the decomposi- 
tion, A is a further parameter of the integration procedure, see section 4. 

2.2.1 A counter-example 

It seems plausible that an arbitrary splitting Xh = Xi + X2 of a Hamiltonian 
vector field need not correspond to a splitting of the Hamiltonian H = H1+H2, 
such that Xi — Xh. for i = 1, 2. Hence the decomposition Xh — X1+X2 does 
not necessarily lead to symplectic integrators. Nevertheless, we will illustrate 
this by an example which is connected with a numerical integrator used for 
bi-partitc spin systems, see [6,9,10]. Such spin systems can be divided into two 
disjoint subsets of spins A and B, such that the interaction is only non-zero 
between spins of different subsets. The first step of the numerical procedure 
consists of fixing the A-spins and calculating the time evolution of all S-spins. 
In the second step the role of A and B is interchanged, and so on. In a single 
step each spin of one subset rotates about the fixed (weighted) sum of all its 
neighboring spins; hence the numerical integrator preserves the volume of the 
total phase space. But, as we will show, this integrator is not symplectic, see 
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also the corresponding remark in [6]. 

It suffices to consider just two spins and a single step of the described numerical 
integrator which solves the equations of motion 

j/i = s.xs„ p, = 0, (17) 



defining a vector field X on V2- We adopt canonical coordinates (pi, Zi,(p2, Z2 
defined in (2) and use the local expression u = dipi A dzi + dip2 A dz2 of the 
symplectic form. After some elementary calculations we obtain 



ixOJ — ipidzi — Zidifi (18) 
= (z2 - ^ly^izif cos((/7i - (/?2)) dzi 

-yj{l-zl){l-zl)sm{ipi - if2)difi . (19) 

Obviously, ix^J is not closed, and hence X does not generate a symplectic 
flow, cf. the discussion after (11). 



3 ;B-partitioned spin systems 



The symplectic integrators considered in section 2.2 are based on a splitting 
of the spin Hamiltonian into a sum of completely integrable Hamiltonians: 
H = J2i Hi. For Heisenberg Hamiltonians 

H{s)^J2 "^M^^^M ■ ' where J^^ e R , (20) 

IJ,<u 



such a splitting is always possible; in fact, each summand in (20) is a com- 
pletely integrable dimcr Hamiltonian. However, it seems favorable to work 
with as few summands as possible, or, equivalently, to work with "large" in- 
tegrable Hamiltonians. To this end we will define a special class of completely 
integrable spin systems called i3-partitioncd systems, following [12]. 
As an example, consider the Heisenberg Hamiltonian of the spin square 

H\J^ Si- S2 + S2- S3 + S3 - S4 + S4- Si . (21) 

It is integrable because it can be written as 

Hn = i ((si + S2 + ss + s^f - (si + ssf - {S2 + s^f) . (22) 
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The grouping of the spins in (22) can be encoded in a "partition tree" 

Ba = {{1, 2, 3, 4}, {1, 3}, {2, 4}, {1}, {2}, {3} {4}} . (23) 



Generahzing this example, we define 

Definition 1 A partition tree B over a finite set {1, . . . , N} is a set of subsets 
of{l,...,N} satisfying 

(1) dl^B and {!,..., N} e B, 

(2) for all M,M' eB either M n M' ^ or M (Z M' or M' G M , 

(3) for all M e B with |M| > 1 there exist Mi,M2 E B such that M = 
M1UM2. 

It follows from definition 1 (2) that the subsets Mi, M2 satisfying M — M1UM2 
in definition 1 (3) are unique, up to their order. Mi, M2 are hence defined for 
all M G jB with |M| > 1. Mi and M2 denote the two uniquely determined 
"branches" starting from M. It follows that i3 is a binary tree with the root 
{1, . . . , A^} and singletons {11} as leaves. More general partitions into k disjoint 
subsets can be reduced to subsequent binary partitions and hence need not be 
considered. For all M e B there is a unique path 

Vm{B) = {M' e B \ M c M'} (24) 



joining M with the root of B. It is linearly ordered since M C M' and M C M" 
imply M' C M" or M" C M' by definition 1 (2). Especially, every element 
// e {!,..., N} belongs to a unique, linearly ordered construction path 

Vi,{B) = {MeB\neM} . (25) 



For fi ^ u E {1, . . . , N} let M^^, G B denote the smallest set of B such that 
/i, z/ G Mfj^^, i. e. M^i, G is the set where both construction paths of fi and ly 
meet the first time. For M 7^ {!,..., N} we will denote by M the "successor" 
of M, that is, the smallest element of Vm{B) except M itself. 

Consider real functions J defined on a partition tree 

J:B — >M. (26) 

satisfying J({/i}) = for all = 1, . . . , A/". Then 

H^Y. JiM^'^K ■ (27) 

IJ,<u 
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defines a Heiscnbcrg Hamiltonian. The corresponding spin system will be 
called a B- partitioned system or sometimes, more precisely, a {B, J) — system. 
For example, the spin square (21) is obtained by the partition tree (23) and 
by the function J with J({1, 2, 3, 4}) = 1 and J(M) = else. 



Let Sm denote the total spin vector of the subsystem M C {1, . . . , N} with 
length Sm- Further, let I]){u,t) denote the 3-dimensional rotation matrix with 
axis cD and angle \uj\ t. In the special case u = 0, Bi{uj,t) denotes the identity 
matrix I. Then the following can be proven, see [12]: 

Theorem 2 Let H he the Hamiltonian of a {B, J) -system. Then its time evo- 
lution is given by 



where the arrow above the product symhol denotes a product according to a 
decreasing sequence of sets M e V^iB) from left to right and J{M) = for 
M = {1,...,N}. 

We note that the time evolution in the presence of a Zeeman term in a Hamil- 

— * — * — * 

tonian of the form H -\- B ■ S, where B is the dimensionless magnetic field, is 
obtained by multiplying (28) from the left with D(S,i). 



If we stick to the symplectic integrators of table 1 and to ^^-partitioned systems 
as completely integrable spin systems, our method only applies to those spin 
systems whose Hamiltonian can be written as the sum of two Hamiltonians of 
B-partitioned subsystems. The spin cube with one additional space diagonal 
is an example which can only be decomposed into at least three S-partitioned 
subsystems. But our method can, in principle, be extended to decompositions 
of the Hamiltonian into more than two summands. As an non-trivial example 
where our method works without modification we mention the spin system of 
= 30 spins which are uniformly coupled according to the edges of an icosi- 
dodecahedron, see [21]. Such a spin system has been physically realized as an 
organic molecule containing 30 paramagnetic Fe-ions, see [22]. In figure 1 the 
planar graph of the icosidodecahedron is decomposed into two i3-partitioned 
subsystems A and B. A consists of 6 disjoint "bow ties" of the form IXI and B 
of 8 disjoint triangles together with 6 single spins. 
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Fig. 1. Decomposition of the graph of the icosidodecahedron into 6 bow ties (dashed 
hnes) and 8 triangles (solid lines). This is the basis of the symplectic integrators 
approximately solving the equations of motion for the corresponding classical spin 
system. 

4 Results 



We have implemented the various symplectic integrators described above us- 
ing the computer algebra software MATHEMATICA 4.0 and have applied 
them to small spin systems. This seems to be sufficient in order to test gen- 
eral properties of the algorithms and to compare the different decompositions 
according to table 1. For more extensive tests and "real life" applications an 
implementation using other computer languages would be advisable. 
For a non-integrable spin system it is impossible to compare the results of 
a numerical integration with the exact result since the exact result is not 
known by definition. Possible tests are observations of conserved quantities 
as the total energy H for non-integrable spin systems or observations of non- 
conserved quantities for integrable spin systems. These tests will be reported 
and discussed in the next subsections. 
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4-1 Total energy 



Figure 2 and 3 show results of numerical integrations of the Hamiltonian equa- 
tions of motion for the spin system corresponding to the icosidodecahedron, see 
figure 1. We choose physical units such that the coupling constant J assumes 
the value 1. For all integrations the time interval is chosen as [0, 100] and the 
time step is A = 0.1. The initial spin configuration is chosen randomly. The 
symplectic integrators applied to this problem are based on a decomposition 
of the icosidodecahedron into 6 bow ties and 8 triangles as explained above. 
Figure 2 shows the total energy and the three components of the total spin 
S as a function of time calculated by the first order integrator STl. Whereas 
S is exactly conserved by all symplectic integrators considered in this article, 
the total energy fluctuates about its initial value with a maximal deviation of 
approximately 7%. 

For symplectic integrators of 4th order the same behavior of the total energy 
can be observed, except that the range of the fluctuation is much smaller. The 
absolute maximal deviation is about 5 ■ 10""^ for FR and 5 ■ 10~^ for OFR and 
ST4, see figure 3. In contrast to these results, a 4th order Runge-Kutta method 
(RK4) yields a systematic drift of the total energy which reaches a deviation of 
1.5 ■ 10-3 at t = 100. This is typical for non-symplectic integrators, see [2], and 
one of the main reasons to adopt symplectic methods for Hamiltonian systems. 
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Fig. 3. Total energy of the icosido decahedron as a function of time calculated by 
ST4, FR, OFR and RK4. 




Fig. 4. Three integrable spin systems used for tests of numerical integrators: The bow 
tie, Nicholas' house and the pyramid. The decomposition into integrable subsystems 
used for symplectic integrators is indicated by solid and dashed lines. 

4-2 Comparison with exact solutions 

We compare non-conserved quantities calculated by the various numerical 
methods with the exact solutions for two integrable systems, the bow tie and 
"Nicholas' house" , see figure 4. The latter is named after a German nursery- 
rhyme ("Das ist das Haus vom Nikolaus"). The time interval [0, 100], the time 
step A = 0.1 and the random choice of the initial configuration is similar as 
in the previous sections. Although the results of the comparison with exact 
solutions shed some light on the respective merits of the different methods, 
it seems dangerous to generalize them to non-integrable problems where the 
distance between near-by solutions may increase exponentially. 
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20 40 60 80 100 

t 

Fig. 5. si • s 5 + S2 • 5*5 + S3 ■ S5 as an exact function of time for the bow tie. 
4.2.1 Bow tie 

Figure 5 shows the quantity s i ■ S5 + S2 • S5 + ss • S5 as an exact function 
of time. Figure 6 shows the absolute deviation 5{si • S5 + S2 • S5 + S3 • S5) 
between the numerical and the exact value in logarithmic scale for the 4th 
order integrators considered above. These deviations seem to increase hnearly 
in time (note the logarithmic scale) but with different orders of magnitude. The 
sharp minima of the logarithmic deviations in this and the following figures 
are due to intersections between the exact and the approximate functions. 
At t — 100 the four integrators can be ordered into a decreasing sequence 
according to their deviations, namely FR, RK4, OFR, ST4, where the ratio 
between two neighbors of this sequence is approximately a factor of 10. It 
is somewhat surprising that the non-symplectic RK4 is better than FR, but 
w r. t. conserved quantities FR should outperform RK4, as shown in section 
4.1. 



4.2.2 Nicholas' house 

It is advisable to consider another example in order to see whether the above 
findings for the bow tie are typical. Figure 7 shows the quantity Si-S3 + S2'S4 for 
the spin system called "Nicholas' house" as an exact function of time. Figure 
8 shows the absolute deviation 6{si ■ S3 + S2 ■ S4) between the numerical and 
the exact value in logarithmic scale for the 4th order integrators considered 
above. These deviations seem to increase again linearly in time (note the 
logarithmic scale). At t = 100 we have two groups, (FR, RK4) and (OFR, 
ST4) with comparable deviations within these groups, where the deviations 
of the second group are almost two orders of magnitude smaller that those of 
the first group. 
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11 RK4 

20 40 60 80 100 

t 

Fig. 6. logxo '^(^1 ■ ■5*5 + •5*2 • s 5 + 5*3 • S5) as a function of time for the bow tie and the 
integrators ST4, FR, OFR, and RK4. 




Fig. 7. s 1 • s 3 + s 2 • S4 as an exact function of time for Nicholas' house. 
4-3 Comparison for given runtime 

From a practical point of view it is not important which numerical proce- 
dure shows the smallest deviations for a fixed time step A but rather for a 
fixed runtime. We will provide a first test of this kind. For this test we have 
to exclude the Runge-Kutta procedures since they are implemented in the 
NDSolve-command of MATHEMATICA and hence their runtime cannot be 
compared with the symplectic integrators programmed in MATHEMATICA 
code. The NDSolve-command of MATHEMATICA 5.0 also allows the choice 
of symplectic integrators, but these integrators are not suited for spin systems 
since they rest on a splitting of the form H = T(p) -|- V^(q), where p, q are 
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Fig. 8. log]^o '^(^1 ■ ^3 + ^2 ■ S4) as a function of time for Nicholas' house and the 
integrators ST4, FR, OFR, and RK4. 




D 20 40 60 80 100 

t 



Fig. 9. s 1 • S5 + s 2 • 5*5 as an exact function of time for the pyramid, 
sets of canonical coordinates. 

The runtime will be measured in terms of the number of "basic operations". 
A basic operation is the calculation of the exact time evolution J-'A{Hi) for 
the Hamiltonian Hi of an integrable subsystem. All basic operations approxi- 
mately require the same cpu-time. The common task is to calculate the quan- 
tity S1-S5 + S2-S5 as a function of t G [0, 100] for a random start configuration of 
the spin pyramid, see figure 4, using maximal 10.000 basic operations. For ev- 
ery numerical procedure the appropriate step size A is separately chosen. The 
results are compared with the exact solution, see figure 9, and the deviations 
are plotted as functions oft in logarithmic scale, see figure 10. The deviations 
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20 40 50 80 100 

t 

Fig. 10. logiQS{si • 5*5 + 5*2 • S5) as a function of time for the pyramid and the 
integrators STl, ST2, ST4, FR, and OFR. 

vary over 8 orders of magnitude and seem to increase with t. It turns out that 
ST4 is three orders of magnitude more precise than ST2 and even five orders 
of magnitude more precise than STl. Whereas FR hes between ST2 and ST4, 
OFR is close to ST4, although its maximal deviation is about two times larger 
than that of ST4. These results indicate that it might be worth while to adopt 
symplectic integrators of even higher order, say, for example, ST6 or ST8. 



5 Summary and outlook 

We have proposed a symplectic integrator scheme for classical spin systems 
based on a splitting of the spin Hamiltonian into two completely integrable 
components corresponding to ^-partitioned subsystems. Further, we have im- 
plemented several variants of this integrator for a selection of small spin sys- 
tems and performed certain tests and comparisons. The results largely conform 
with the expectations; an interesting finding is that, for fixed runtime, higher 
order algorithms yield marked improvements of the precision. This accords 
with the results of [2], section V.3.2, where, however, no further improvement 
occurs beyond the order of 8. 

Of course, these tests are only preliminary and should be extended to include, 
for instance, more spin systems, the longtime behavior and the influence of 
different decompositions of the Hamiltonian. Also we have not compared our 
method with other methods which are energy- and volume-preserving, but not 
symplectic [6,9,10]. Our method cannot be applied to an arbitrary Hamilto- 
nian spin system without taking additional measures. This is a draw-back, but 
simultaneously an advantage since it means that one has to adapt the method 
for a given system in order to find an optimal algorithm. In view of the ap- 
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plicability the perhaps most pressing generahzation would be to consider the 
case of more than two integrable components of the Hamiltonian. 
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